## FUNCTIONS

mean.IBD <- function (data){
#	res <- as.numeric(data[4])+as.numeric(data[5])*0.5 # mean IBD = p(IBD=2) + p(IBD=1) * 0.5
	res <- as.numeric(data[2])+as.numeric(data[3])*0.5 # mean IBD = p(IBD=2) + p(IBD=1) * 0.5
}

abs.diff <- function(data){
#	res <- abs(as.numeric(data[7])*2 - as.numeric(data[8]))
	res <- abs(as.numeric(data[5])*2 - as.numeric(data[6]))
}

ratio <- function(data){
	res <- log((as.numeric(data[5])*2)/as.numeric(data[6]))
}

#######


# arguments
args <- commandArgs(TRUE)

kinship <- read.table(file=args[1])
#mean pair kinship for all chr
meanPairKinship <- tapply(X=kinship$V6,INDEX=list(kinship$V2),FUN=mean) 
pairs <- as.character(names(meanPairKinship))
meanpkinship <- cbind(pairs,as.numeric(meanPairKinship))

# mean individual kinship for all chr
meanIndivKinship <- tapply(X=kinship$V6,INDEX=list(kinship$V3),FUN=mean)
indiv <- as.character(names(meanIndivKinship))
meanikinship <- cbind(indiv,as.numeric(meanIndivKinship))

# mean pair IBD
reskin1 <- read.table(file=args[2])
tmpreskin1 <- reskin1[-c(1,2,3)]
names(tmpreskin1) <- c("V1","V2","V3","V4")
mergep1 <- merge(tmpreskin1,meanpkinship,by.x="V1",by.y="pairs")

reskin2 <- read.table(file=args[3])
tmpreskin2 <- reskin2[-c(1,2,3)]
names(tmpreskin2) <- c("V1","V2","V3","V4")
mergep2 <- merge(tmpreskin2,meanpkinship,by.x="V1",by.y="pairs")

mergep <- rbind(mergep1,mergep2)
tablep <- unique(mergep)

# difference : mean kinship*2 - mean IBD (observed - expected)
tablep$meanIBD <- apply(X=tablep,MARGIN=1,FUN=mean.IBD)
#tmptablep <- tablep[-c(2,3)]
tablep$diff <- apply(X=tablep,MARGIN=1,FUN=abs.diff)
names(tablep) <- c("pairs","pIBD2","pIBD1","pIBD0","kinshipObs","meanIBDexpected","abs_diff_kinshipIBD")

# mean individuals IBD
reskin <- rbind(reskin1,reskin2)
tmpreskin <- reskin[-c(1,3,4)]
tmpreskin$meanIBD<- apply(X=tmpreskin,MARGIN=1,FUN=mean.IBD)
meanIBD <- tapply(X=tmpreskin$meanIBD,INDEX=list(tmpreskin$V2),FUN=mean)
#meanIBD <- tapply(X=reskin$meanIBD,INDEX=list(reskin$V2),FUN=mean)
indiv <- as.character(names(meanIBD))
meanIndivIBD <- cbind(indiv,as.numeric(meanIBD))
tablei <- merge(meanIndivIBD,meanikinship,by="indiv")
fidiid <- unique(reskin1[c(1,2)])
tmptablei <- merge(tablei,fidiid,by.x="indiv",by.y="V2")
tablei <- tmptablei[c(4,1,2,3)]
tablei$tmp1 <-0
tablei$tmp2 <-0
tablei <- tablei[c(1,2,5,6,4,3)]

# difference
tablei$diff_kinshipIBD <- apply(X=tablei,MARGIN=1,FUN=abs.diff)
tablei$abs_ratio <- apply(X=tablei,MARGIN=1,FUN=ratio) 
tablei <- tablei[c(1,2,5,6,7,8)]
names(tablei) <- c("FID","IID","meanKinship","meanIBD","abs_diff_kinshipIBD","log_diff_ratio")



#indivmeandiff  <- tapply(X=table1$diff,INDEX=list(table1$V3),FUN=mean)
#indivmeandiff <- tapply(X=table$diff,INDEX=table$V3,FUN=mean)
#indiv <- as.character(names(indivmeandiff))
#meandiff <- cbind(indiv,indivmeandiff)
#meandiff <- na.omit(meandiff)

# save datas
write.table(tablep,file=paste(args[4],"_meanPairKinship_meanIBD.txt",sep=""),quote=F,row.names=F)
write.table(tablei,file=paste(args[4],"_meandiff_individuals.txt",sep=""),quote=F,row.names=F)

